| Gen | Muestras | |||||
|---|---|---|---|---|---|---|
| sample1 | sample2 | sample3 | sample4 | sample5 | sample6 | |
| gene1 | 85 | 76 | 103 | 107 | 140 | 124 |
| gene2 | 1 | 6 | 11 | 6 | 1 | 4 |
| gene3 | 80 | 98 | 39 | 82 | 97 | 113 |
| gene4 | 92 | 83 | 59 | 85 | 100 | 98 |
| gene5 | 36 | 24 | 18 | 50 | 22 | 15 |
| gene6 | 0 | 0 | 1 | 4 | 2 | 3 |
1 Análisis de expresión génica diferencial (RNA-Seq)
El análisis de la expresión génica diferencial (RNA-Seq) es una de las técnicas bioinformáticas que más auge ha experimentado en la última década, ya que permite cuantificar la expresión de un gen a nivel celular y comparar estas expresiones entre distintos grupos experimentales para identificar asociaciones entre genes y factores experimentales (como enfermedades o tratamientos).
El proceso de un análisis genético diferencial tiene distintas etapas, que van desde la obtención del las células a estudiar, pasando por la secuenciación del RNA que contienen, hasta el análisis estadístico de estas frecuencias para determinar el nivel de expresión de cada gen y la comparación entre los grupos experimentales.


En este tutorial nos centraremos en el análisis estadístico de las frecuencias de expresión de los genes, que normalmente requiere las siguientes etapas:
Creación de la estructura de datos con las frecuencias de expresión génica.
Preprocesamiento de datos.
Análisis exploratorio de los datos.
Análisis de expresión génica diferencial.
Visualización de resultados.
Existen multitud de paquetes en el repositorio Bioconductor para realizar el análisis de expresión génica diferencial. En este tutorial veremos dos de los más usados: EdgeR y DesSeq2.
1.1 Análisis de expresión génica diferencial con edgeR
El paquete edgeR es uno de los principales paquetes usados para el análisis de expresión génica diferencial que incorpora funciones para todas las etapas del análisis estadístico. Está disponible en el repositorio Bioconductor.
A continuación se explica cómo realizar las distintas etapas del análisis estadístico con este paquete.
1.1.1 Estructura de datos
El análisis estadístico comienza siempre a partir de la matriz con las frecuencias o conteos de expresión de cada gen en el conjunto de muestras analizadas. Las filas de esta matriz corresponden a los genes observados y las columnas a las muestras analizadas.
Ejemplo 1.1 A continuación se muestran las primeras filas de una matriz de frecuencias de expresión génica. Cada casilla recoge el número de veces que el gen de la fila correspondiente se ha observado en la muestra de la columna correspondiente.
El paquete edgeR utiliza la estructura de datos DGEList para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se utiliza la función
DGEList(frec, group = grupo): Crea una estructura del tipoDGEListcon la matriz de frecuencias de expresión génicafrec(con genes en filas y muestras en columnas) y el vectorgrupocon los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.
library(edgeR)
library(DEFormats)
frec <- simulateRnaSeqData()
grupo <- rep(c("A", "B"), each = 3)
dge <- DGEList(frec, group = grupo)
dgeAn object of class "DGEList"
$counts
sample1 sample2 sample3 sample4 sample5 sample6
gene1 85 76 103 107 140 124
gene2 1 6 11 6 1 4
gene3 80 98 39 82 97 113
gene4 92 83 59 85 100 98
gene5 36 24 18 50 22 15
995 more rows ...
$samples
group lib.size norm.factors
sample1 A 42832 1
sample2 A 40511 1
sample3 A 39299 1
sample4 B 43451 1
sample5 B 40613 1
sample6 B 43662 1
Esta estructura de datos es una lista con dos atributos. El atributo counts contiene la matriz de frecuencias de expresión génica, y el atributo samples es un data frame con información sobre las muestras estudiadas. Cada fila de este data frame se corresponde con una columna de la matriz de frecuencias y contiene las siguientes columnas
| Columna | Descripción |
|---|---|
group |
Grupo experimental al que pertenece la muestra. |
lib.size |
tamaño de la librería (suma de frecuencias de la muestra). |
norm.factors |
Factor de normalización. |
La estructura de datos DGEList puede contener opcionalmente el atributo genes con anotaciones de los genes observados en las filas de la matriz de frecuencias.
Ejemplo 1.2 Veamos cómo crear la estructura de datos correspondiente al estudio de Sheridan JM, et al. (2015) donde se obtuvieron datos de tres poblaciones de células (basal, progenitor luminal (LP) y luminal maduro (ML)) seleccionadas de las glándulas mamarias de ratones vírgenes hembra, cada una por triplicado. Los datos pueden obtenerse del repositorio Gene Expression Omnibus (GEO) mediante el identificador GSE63310.
En primer lugar cargamos la matriz de frecuencias.
frecuencias <- read.csv("datos/matriz-frecuencias-genes.csv", row.names = 1)
head(frecuencias) LP1 ML1 B1 B2 ML2 LP2 B3 ML3 LP3
497097 1 2 342 526 3 3 535 2 0
100503874 0 0 5 6 0 0 5 0 0
100038431 0 0 0 0 0 0 1 0 0
19888 0 1 0 0 17 2 0 1 0
20671 1 1 76 40 33 14 98 18 8
27395 431 771 1368 1268 1564 769 818 468 342
Después cargamos los datos de las muestras con los factores experimentales.
muestras <- read.csv("datos/muestras.csv")
muestras Id Grupo Lote
1 LP1 LP L004
2 ML1 ML L004
3 B1 Basal L004
4 B2 Basal L006
5 ML2 ML L006
6 LP2 LP L006
7 B3 Basal L006
8 ML3 ML L008
9 LP3 LP L008
A continuación creamos la estructura de datos DGEList.
# Creación de la estructura de datos DGEList.
dge <- DGEList(counts = frecuencias, group = muestras$Grupo)
# Añadimos también el Lote al dafa frame de las muestras.
dge$samples$Lote <- muestras$Lote
dgeAn object of class "DGEList"
$counts
LP1 ML1 B1 B2 ML2 LP2 B3 ML3 LP3
497097 1 2 342 526 3 3 535 2 0
100503874 0 0 5 6 0 0 5 0 0
100038431 0 0 0 0 0 0 1 0 0
19888 0 1 0 0 17 2 0 1 0
20671 1 1 76 40 33 14 98 18 8
27174 more rows ...
$samples
group lib.size norm.factors Lote
LP1 LP 32863052 1 L004
ML1 ML 35335491 1 L004
B1 Basal 57160817 1 L004
B2 Basal 51368625 1 L006
ML2 ML 75795034 1 L006
LP2 LP 60517657 1 L006
B3 Basal 55086324 1 L006
ML3 ML 21311068 1 L008
LP3 LP 19958838 1 L008
A continuación anotamos los genes de las filas de la matriz de frecuencias. Para ello debe utilizarse un paquete específico para el genoma de la especie de donde proviene el RNA (Mus.muculus en este caso para el genoma del ratón).
library(Mus.musculus)
geneid <- rownames(dge)
genes <- select(Mus.musculus, keys = geneid, columns = c("SYMBOL", "TXCHROM"), keytype = "ENTREZID")
# Eliminamos duplicidad de algunos genes
genes <- genes[!duplicated(genes$ENTREZID),]
dge$genes <- genes
dgeAn object of class "DGEList"
$counts
LP1 ML1 B1 B2 ML2 LP2 B3 ML3 LP3
497097 1 2 342 526 3 3 535 2 0
100503874 0 0 5 6 0 0 5 0 0
100038431 0 0 0 0 0 0 1 0 0
19888 0 1 0 0 17 2 0 1 0
20671 1 1 76 40 33 14 98 18 8
27174 more rows ...
$samples
group lib.size norm.factors Lote
LP1 LP 32863052 1 L004
ML1 ML 35335491 1 L004
B1 Basal 57160817 1 L004
B2 Basal 51368625 1 L006
ML2 ML 75795034 1 L006
LP2 LP 60517657 1 L006
B3 Basal 55086324 1 L006
ML3 ML 21311068 1 L008
LP3 LP 19958838 1 L008
$genes
ENTREZID SYMBOL TXCHROM
1 497097 Xkr4 chr1
2 100503874 Gm19938 <NA>
3 100038431 Gm10568 <NA>
4 19888 Rp1 chr1
5 20671 Sox17 chr1
27174 more rows ...
1.1.2 Preprocesamiento
Una vez preparada la estructura de datos, el siguiente paso es el preprocesamiento de datos. Normalmente comprende las siguientes tareas:
- Normalización de las frecuencias.
- Eliminación de los genes con poca expresión.
- Normalización de las distribuciones de frecuencias de expresión génicas.
1.1.2.1 Normalización de las frecuencias
El objetivo principal es normalizar las frecuencias de expresión génica para eliminar el efecto sobre las frecuencias de factores como la profundidad de secuenciado (a mayor profundidad de secuenciado mayores frecuencias) o el tamaño de los genes (a mayor tamaño mayores frecuencias). Generalmente se usan las siguientes transformaciones
Frecuencias por millón (CPM). Se divide cada frecuencia por el tamaño en millones de la librería de la muestra a la que pertenece. Se realiza con la función
cpm(dge).Logaritmo en base 2 de las frecuencias por millón (log2-CPM). Se calcula a partir de la anterior mediante la fórmula \(\log_2(CPM+\frac{2}{L})\), donde \(L\) es la longitud media de las librerías en millones. El término \(\frac{2}{L}\) que se añade permite calcular el logaritmo para frecuencias nulas (ya que el logaritmo de 0 no existe). Se realiza con la función
cmp(dge, log = TRUE).Lecturas por kilobase de transcripción por millón (RPKM): Se realiza con la función
rpkm(dge, longitud_genes)pasando la longitud de los genes.Fragmentos por kilobase de transcripción por millón (FPKM).
Las dos primeras no tienen en cuenta el tamaño de los genes, pero suelen usarse para la expresión génica diferencial donde el tamaño de los genes se supone constante en las distintas muestras.
Ejemplo 1.3 Siguiendo con el ejemplo anterior, calculamos las frecuencias por millón y el logaritmo de las frecuencias por millón.
cpm <- cpm(dge)
summary(cpm) LP1 ML1 B1 B2
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.578 Median : 0.736 Median : 0.892 Median : 0.895
Mean : 36.793 Mean : 36.793 Mean : 36.793 Mean : 36.793
3rd Qu.: 19.536 3rd Qu.: 23.546 3rd Qu.: 24.221 3rd Qu.: 23.341
Max. :27807.947 Max. :11546.719 Max. :7951.408 Max. :7389.433
ML2 LP2 B3 ML3
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.699 Median : 0.711 Median : 0.908 Median : 0.845
Mean : 36.793 Mean : 36.793 Mean : 36.793 Mean : 36.793
3rd Qu.: 23.827 3rd Qu.: 19.928 3rd Qu.: 21.439 3rd Qu.: 24.260
Max. :7955.680 Max. :29572.361 Max. :9376.973 Max. :7865.350
LP3
Min. : 0.000
1st Qu.: 0.000
Median : 0.752
Mean : 36.793
3rd Qu.: 21.594
Max. :16500.710
lcpm <- cpm(dge, log = TRUE)
summary(lcpm) LP1 ML1 B1 B2
Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
Median :-0.6847 Median :-0.3589 Median :-0.09513 Median :-0.0901
Mean : 0.1714 Mean : 0.3312 Mean : 0.43559 Mean : 0.4089
3rd Qu.: 4.2913 3rd Qu.: 4.5601 3rd Qu.: 4.60081 3rd Qu.: 4.5475
Max. :14.7632 Max. :13.4952 Max. :12.95700 Max. :12.8513
ML2 LP2 B3 ML3
Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
Median :-0.4281 Median :-0.4064 Median :-0.07152 Median :-0.1704
Mean : 0.3225 Mean : 0.2529 Mean : 0.40428 Mean : 0.3708
3rd Qu.: 4.5772 3rd Qu.: 4.3199 3rd Qu.: 4.42513 3rd Qu.: 4.6031
Max. :12.9578 Max. :14.8520 Max. :13.19491 Max. :12.9413
LP3
Min. :-4.5074
1st Qu.:-4.5074
Median :-0.3300
Mean : 0.2749
3rd Qu.: 4.4355
Max. :14.0102
1.1.2.2 Eliminación de genes con poca expresión
Aunque el objetivo del análisis de la expresión génica diferencial es detectar los genes que se expresan en un grupo experimental en comparación con los otros, cuando un gen no se expresa en ninguna de las muestras no tiene interés para el estudio y puede eliminarse.
Existen diferentes criterios para eliminar los genes con poca expresión. El paquete EdgeR incorpora la siguiente función
filterByExpr(dge): Devuelve un vector de booleanos con nombres correspondientes a los genes de la estructura DEGListdge. Por defecto devuelveTRUEpara cualquier gen con una frecuencia mayor o igual que 10 en al menos un número de muestras igual al del menor grupo experimental.
Ejemplo 1.4 Veamos cuántos genes no tienen expresión en ninguna muestra del ejemplo anterior.
# Número de genes con expresión nula.
table(rowSums(dge$counts) == 0)
FALSE TRUE
22026 5153
El siguiente gráfico muestra la distribución del logaritmo en base 2 de las frecuencias por millón.
# Definimos una función para dibujar la distribución del logaritmo de las frecuencias por millón.
dist_logcpm <- function(lcpm) {
lcpm |>
as.tibble() |>
pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
ggplot(aes(x = LogCPM, color = Muestra)) +
geom_density() +
labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}
dist_logcpm(lcpm)Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Como se observa ha una porcentaje significativo de genes que tiene una expresión muy baja (valores negativos).
A continuación eliminamos de la estructura de datos los genes con poca expresión.
filtro <- filterByExpr(dge)
dge <- dge[filtro, , keep.lib.sizes = F]Y volvemos a dibujar la distribución del logaritmo en base 2 de las frecuencias por millón.
lcpm <- cpm(dge, log = TRUE)
dist_logcpm(lcpm)
Como se aprecia en el diagrama, el porcentaje de genes con baja expresión ha disminuido significativamente.
1.1.2.3 Normalización de las distribuciones de frecuencias de expresión génicas
Durante el proceso de secuenciación puede haber factores externos no biológicos que afecten a la expresión de los genes. Por ejemplo, la muestras del primer lote pueden tener mayores frecuencias que las del segundo lote. Como se supone que las muestras replicadas deben tener un rango y una distribución de frecuencias similares, en esta etapa se aplica otro procedimiento de normalización para garantizar que la distribución de frecuencias de cada muestra es similar. Para ello el paquete edgeR utiliza el método de las medias recortadas de los valores M por medio de la función
calcNormFactors(dge, method = "TMM"): Calcular los factores de escalado del tamaño de las librerías de cada muestra. Estos factores se guardan automáticamente en el data frame con la información de las muestrasdge$samples.
Ejemplo 1.5 Siguiendo con el mismo ejemplo, calculamos los factores de escalado para cada muestra.
dge <- calcNormFactors(dge, method = "TMM")
dge$samples group lib.size norm.factors Lote
LP1 LP 32857304 0.8943956 L004
ML1 ML 35328624 1.0250186 L004
B1 Basal 57147943 1.0459005 L004
B2 Basal 51356800 1.0458455 L006
ML2 ML 75782871 1.0162707 L006
LP2 LP 60506774 0.9217132 L006
B3 Basal 55073018 0.9961959 L006
ML3 ML 21305254 1.0861026 L008
LP3 LP 19955335 0.9839203 L008
A continuación se muestran los diagramas de cajas de las distribuciones normalizadas tras aplicar los factores de escalado.
box_logcpm <- function(dge){
muestras <- dge$samples
muestras$Muestra <- row.names(muestras)
cpm(dge, log = TRUE) |>
as.tibble() |>
pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
left_join(muestras, by = "Muestra") |>
ggplot(aes(x = Muestra, y = LogCPM, fill = group)) +
geom_boxplot() +
labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}
box_logcpm(dge)
Como se puede apreciar, después de la normalización, todas las muestras presentan una distribución de frecuencias similar.
1.1.3 Análisis exploratorio
Una vez preprocesados los datos comienza el análisis estadístico propiamente dicho. En una primera fase se suele realizar un análisis exploratorio de los datos que suele incluir los siguientes análisis:
- Escalado multidimensional (análisis de componentes principales).
1.1.3.1 Escalado multidimensional
El escalado multidimensional mediante componente principales permite ver qué muestras son similares en cuando a la distribución de frecuencias de expresión génica. Los componentes principales son combinaciones lineales de los genes de la matriz de frecuencias que son ortogonales entre sí. El primer componente principal recoge la dimensión con mayor variabilidad de las frecuencias. El segundo recoge la dimensión don la mayor variabilidad no explicada por el primer componente principal y así sucesivamente. Normalmente los dos primeros componentes principales suelen recoger un porcentaje bastante alto de la variabilidad de las frecuencias. Al representar las muestras en estas dos dimensiones, las muestras más próximas entre sí, presentan una distribución de frecuencias de expresión génica similar. Para realizar esta representación se puede utilizar la siguiente función del paquete limma:
plotMDS(dge): Realiza un diagrama de componentes principales de la matriz de frecuencias de la estructura de datosdge.
Alternativamente, se pueden calcular los componentes principales mediante la función prcomp del paquete stats y luego usar la función autoplot del paquete ggfortify para dibujar el diagrama de componentes principales.
Ejemplo 1.6 A continuación se muestra cómo obtener el diagrama de componentes principales de la matriz de frecuencias de nuestro ejemplo.
plotMDS(dge, col = as.numeric(dge$samples$group), main = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")
library(ggfortify)
autoplot(prcomp(t(lcpm)), data = dge$samples, color = "group", loadings = F, loadings.label = F) +
labs(title = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")
Como se puede apreciar las muestras de cada grupo experimental aparecen agrupadas. La mayor diferencia (a lo largo del primer componente principal) se da entre el grupo basal y los otros dos grupos. Esto indica que cuando se realice el análisis de expresión génica diferencial habrá bastantes genes con diferencias de expresión significativa entre el grupo basal y los otros dos, mientras que no habrá tantos al comparar los grupos LP y ML.
Otra opción muy interesante es el paquete Glimma que permite dibujar un diagrama de componentes principales interactivo mediante la función
glMDSPlot(lcpm): Dibuja un diagrama de componentes principales interactivo de la matriz de frecuenciaslcpm.
Ejemplo 1.7 A continuación se muestra cómo obtener el diagrama de componentes principales para nuestro ejemplo con el paquete Glimma.
library(Glimma)
library(DESeq2)
#glMDSPlot(lcpm, groups = dge$samples[,c(2,5)])
dds <- DESeqDataSetFromMatrix(
countData = dge$counts,
colData = dge$samples,
rowData = dge$genes,
design = ~group
)
glimmaMDS(dds)1.1.4 Análisis de expresión génica diferencial
La última etapa del análisis, y la más importante, es la detección de los genes con una expresión significativamente diferente entre los grupos experimentales. Para ello se suelen utilizar modelos lineales o modelos lineales generalizados. En general, la estimación de los parámetros del modelo ajustado depende de la distribución teórica usada para modelizar las frecuencias de expresión génica. El paquete limma, por ejemplo, realiza un ajuste de modelo lineal suponiendo que las distribuciones de las variables implicadas son normales, mientras que el paquete edgeR modeliza las frecuencias de expresión de los genes observadas mediante la distribución binomial negativa, que es mucho más realista.
En general en las distribuciones de frecuencias de expresión génica, se ha observado que la varianza no es independiente de la media. Los métodos que modelizan las frecuencias mediante el modelo binomial negativo asumen una relación cuadrática entre la varianza y media.
El siguiente paso es estimar la dispersión de las frecuencias de expresión génica para cada gen. Para ello edgeR utiliza el método de la máxima verosimilitud condicionada y ajustada por percentiles, mediante la función
estimateDisp(dge, diseño): Realiza una estimación de la dispersión de las frecuencias de expresión génica contenidas en la estructura del tipo DGEList dad endge, de acuerdo al diseño experimental dado endiseño.
Ejemplo 1.8
dge <- estimateDisp(dge)Using classic mode.
Una vez ajustado el modelo Binomial Negativo y estimada la dispersión, solo queda aplicar el contraste de comparación de la expresión génica diferencial entre los grupos experimentales. Para ello edgeR utiliza el test exacto similar al test exacto de Fisher, pero adaptado a la distribución Binomial Negativa. Para ello se usa la función
exacTest(dge, pair=grupos): Realiza el contraste de comparación de la expresión génica a partir del las frecuencias de aparción de genes contenidas en la estructura del tipo DGEList dada endge, entre los grupos experimentales dados en el vectorgrupos(solo puede contener dos grupos).
Ejemplo 1.9 En nuestro ejemplo realizamos tres contrastes para todos los posibles pares de grupos experimentales.
exact_B_LP <- exactTest(dge, pair= c("Basal", "LP"))
exact_B_ML <- exactTest(dge, pair= c("Basal", "ML"))
exact_LP_ML <- exactTest(dge, pair= c("LP", "ML"))Tras realizar el contraste se puede utilizar la siguiente función para obtener un resumen con el número de genes subrexpresados y sobrexpresados significativamente.
summary(decideTests(test)): Devuelve una tabla con los genes subexpresados, sobreexpresados significativamente, así como lo que no presentan cambios significativos a partir de los resultados del contrastetest.
Finalmente, para ver los genes con diferencias de expresión más estadísticamente significativas se puede utilizar la función
topTags(test): Devuelve un data frame con los genes con diferencias de expresión estadísticamente más significativas entre los grupos experimentales comparados en objeto de resultados del contrastetest.
Ejemplo 1.10 Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo Basal y LP.
summary(decideTests(exact_B_LP)) LP-Basal
Down 4864
NotSig 7203
Up 4557
Genes con diferencias más significativas entre el grupo Basal y el grupo LP.
as.data.frame(topTags(exact_B_LP)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | PValue | FDR |
|---|---|---|---|---|---|---|---|
| 67451 | 67451 | Pkp2 | chr16 | 5.661601 | 5.672839 | 2.184632e-164 | 3.631733e-160 |
| 19253 | 19253 | Ptpn18 | chr1 | 5.588760 | 5.344834 | 3.519102e-160 | 2.925078e-156 |
| 218518 | 218518 | Marveld2 | chr13 | 5.216899 | 6.054980 | 4.369286e-157 | 2.421167e-153 |
| 50722 | 50722 | Dkkl1 | chr7 | 6.296222 | 6.098528 | 6.043793e-156 | 2.511800e-152 |
| 227960 | 227960 | Gca | chr2 | 5.516041 | 5.060151 | 6.024796e-144 | 2.003124e-140 |
| 242505 | 242505 | Rasef | chr4 | 6.071991 | 6.704767 | 2.430700e-143 | 6.734659e-140 |
| 320007 | 320007 | Sidt1 | chr16 | 6.297688 | 4.897542 | 1.142473e-135 | 2.713210e-132 |
| 22249 | 22249 | Unc13b | chr4 | 4.343846 | 6.600123 | 1.313443e-134 | 2.729335e-131 |
| 66871 | 66871 | Cpne8 | chr15 | -4.685486 | 5.862087 | 4.580752e-129 | 8.461158e-126 |
| 269132 | 269132 | Colgalt2 | chr1 | -4.989659 | 5.011534 | 4.878508e-126 | 8.110031e-123 |
Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo Basal y ML.
summary(decideTests(exact_B_ML)) ML-Basal
Down 4773
NotSig 7011
Up 4840
Genes con diferencias más significativas entre el grupo Basal y el grupo ML.
as.data.frame(topTags(exact_B_ML)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | PValue | FDR |
|---|---|---|---|---|---|---|---|
| 50722 | 50722 | Dkkl1 | chr7 | 6.728394 | 6.098528 | 4.166088e-173 | 6.925705e-169 |
| 242505 | 242505 | Rasef | chr4 | 6.678297 | 6.704767 | 2.371625e-165 | 1.971295e-161 |
| 21804 | 21804 | Tgfb1i1 | chr7 | -6.120166 | 6.012779 | 9.151423e-156 | 5.071108e-152 |
| 22249 | 22249 | Unc13b | chr4 | 4.561206 | 6.600123 | 3.599522e-150 | 1.495961e-146 |
| 20661 | 20661 | Sort1 | chr3 | 4.948551 | 7.722875 | 8.615672e-144 | 2.864539e-140 |
| 12521 | 12521 | Cd82 | chr2 | 4.666737 | 8.007651 | 8.771904e-143 | 2.430402e-139 |
| 67451 | 67451 | Pkp2 | chr16 | 5.100406 | 5.672839 | 7.763029e-140 | 1.843608e-136 |
| 218518 | 218518 | Marveld2 | chr13 | 4.815118 | 6.054980 | 7.247449e-139 | 1.506020e-135 |
| 78896 | 78896 | Ecrg4 | chr1 | -6.587115 | 5.670792 | 9.526597e-139 | 1.759668e-135 |
| 214968 | 214968 | Sema6d | chr2 | -5.935614 | 6.101960 | 1.913875e-125 | 3.181626e-122 |
Resumen con los genes subexpresados y sobreexpresados significativamente entre el grupo LP y ML.
summary(decideTests(exact_LP_ML)) ML-LP
Down 2432
NotSig 11226
Up 2966
Genes con diferencias más significativas entre el grupo LP y el grupo ML.
as.data.frame(topTags(exact_LP_ML)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | PValue | FDR |
|---|---|---|---|---|---|---|---|
| 94212 | 94212 | Pag1 | chr3 | -2.616277 | 5.916572 | 9.314994e-31 | 1.548525e-26 |
| 214968 | 214968 | Sema6d | chr2 | -2.378756 | 6.101960 | 3.482661e-27 | 2.894788e-23 |
| 114479 | 114479 | Slc5a5 | chr8 | 2.879278 | 7.770145 | 1.246288e-26 | 6.906100e-23 |
| 207592 | 207592 | Tbc1d16 | chr11 | 2.308784 | 5.659683 | 6.219699e-26 | 2.584907e-22 |
| 13617 | 13617 | Ednra | chr8 | -3.187364 | 4.014091 | 1.557011e-25 | 5.176750e-22 |
| 231605 | 231605 | Galnt9 | chr5 | 3.289480 | 2.482287 | 2.766785e-24 | 7.665838e-21 |
| 12614 | 12614 | Celsr1 | chr15 | 2.826770 | 5.693953 | 1.525133e-23 | 3.621973e-20 |
| 14778 | 14778 | Gpx3 | chr11 | 3.045448 | 9.996423 | 4.769484e-23 | 9.910988e-20 |
| 107895 | 107895 | Mgat5 | chr1 | 2.025411 | 5.535760 | 9.509983e-23 | 1.756600e-19 |
| 320127 | 320127 | Dgki | chr6 | 2.415427 | 4.951953 | 1.261948e-22 | 2.097863e-19 |
Otra alternativa, para diseños experimentales que incluye más de un factor, es usar modelos lineales generalizados (GLM). El primer paso es definir la matriz del diseño del experimento, que incluye las variables que definen grupos experimentales. Para ello se utiliza la función
model.matrix(formula-diseño): Construye una matriz con el diseño experimental de acuerdo a la fórmula dada enformula-diseño. Esta fórmula es similar a la fórmula que que utiliza en un ANOVA para definir el diseño experimental, donde las variables se añaden con el operador+, y se hacen interactuar con el operador*.
Ejemplo 1.11 Continuando con el mismo ejemplo, definimos la matriz de diseño del modelo.
diseño <- model.matrix(~ 0 + dge$samples$group + dge$samples$Lote)
colnames(diseño) <- gsub("dge\\$samples\\$group", "", colnames(diseño))
colnames(diseño) <- gsub("dge\\$samples\\$", "", colnames(diseño))
diseño Basal LP ML LoteL006 LoteL008
1 0 1 0 0 0
2 0 0 1 0 0
3 1 0 0 0 0
4 1 0 0 1 0
5 0 0 1 1 0
6 0 1 0 1 0
7 1 0 0 1 0
8 0 0 1 0 1
9 0 1 0 0 1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"
attr(,"contrasts")$`dge$samples$Lote`
[1] "contr.treatment"
Al igual que antes, antes de realizar el contraste hay que estimar la dispersión conjunta y para cada gen, pero ahora hay que introducir el diseño como un parámetro de la función estimateDisp.
Ejemplo 1.12 Realizamos la estimación de la dispersión para nuestro ejemplo.
dge <- estimateDisp(dge, diseño)Una vez estimada la dispersión y ajustados los modelos generalizados binomiales negativos, ya se puede realizar el contraste de comparación de expresión génica. Para ello se utilizan las siguientes funciones
glmQLFit(dge, diseño): Realiza el ajuste del modelo generalizado de comparación de expresión génica mediante el test F de cuasi verosimilitud, con las frecuencias de expresión de la estructura del tipo DGEList dada endgey de acuerdo al diseño experimental indicado endiseño.makeContrast(grupo1 - grupo2, levels=diseño): Define los gruposgrupo1ygrupo2a comparar en el contraste pares para la diferencia en la expresión génica.glmQLFTest(modelo-ajustado, contrast=contraste). Realiza el contraste de comparación por pares a partir del modelo ajustadomodelo-ajustado. El modelo ajustado tiene varios coeficientes dependiendo del número grupos experimentales, el primero corresponde al ajuste base para el grupo control, y los sucesivos a las diferencias de los otros grupos experimentales con el control.glmTreat(modelo-ajustado, contrast=contraste, lfc = n). Realiza el contraste de comparación por pares similar al de la función anterior, pero descartando los genes con un logaritmo en base dos de la razón de cambio (logFC) menor que el valor dado en el parámetrolfc.
Ejemplo 1.13 Realizamos el ajuste del modelo linear generalizado para nuestro ejemplo.
modelo <- glmQLFit(dge, diseño)Y, a continuación, los contrastes por pares. En primer lugar, comparamos Basal con LP y mostramos los genes más diferenciados.
contraste <- makeContrasts(Basal - LP, levels = diseño)
glmQL_B_LP <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_B_LP)) 1*Basal -1*LP
Down 4523
NotSig 7192
Up 4909
as.data.frame(topTags(glmQL_B_LP)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | F | PValue | FDR |
|---|---|---|---|---|---|---|---|---|
| 242505 | 242505 | Rasef | chr4 | -5.941271 | 6.704195 | 1863.844 | 5.629045e-11 | 4.877025e-07 |
| 67451 | 67451 | Pkp2 | chr16 | -5.745868 | 5.671932 | 1845.036 | 5.867451e-11 | 4.877025e-07 |
| 19253 | 19253 | Ptpn18 | chr1 | -5.655466 | 5.343492 | 1616.133 | 1.008466e-10 | 5.588246e-07 |
| 53624 | 53624 | Cldn7 | chr11 | -5.527387 | 7.529234 | 1327.691 | 2.251566e-10 | 6.258841e-07 |
| 14275 | 14275 | Folr1 | chr7 | -6.925474 | 5.510509 | 1313.864 | 2.349910e-10 | 6.258841e-07 |
| 218518 | 218518 | Marveld2 | chr13 | -5.153654 | 6.054055 | 1300.770 | 2.448002e-10 | 6.258841e-07 |
| 70350 | 70350 | Basp1 | chr15 | -6.086771 | 6.624263 | 1223.255 | 3.145873e-10 | 6.258841e-07 |
| 228543 | 228543 | Rhov | chr2 | -6.263337 | 6.939962 | 1112.552 | 4.632731e-10 | 6.258841e-07 |
| 70737 | 70737 | Cgn | chr3 | -5.466154 | 6.644032 | 1075.182 | 5.325457e-10 | 6.258841e-07 |
| 320007 | 320007 | Sidt1 | chr16 | -6.345681 | 4.895659 | 1071.587 | 5.398676e-10 | 6.258841e-07 |
A continuación comparamos Basal con ML.
contraste <- makeContrasts(Basal - ML, levels = diseño)
glmQL_B_ML <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_B_ML)) 1*Basal -1*ML
Down 4811
NotSig 7060
Up 4753
as.data.frame(topTags(glmQL_B_ML)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | F | PValue | FDR |
|---|---|---|---|---|---|---|---|---|
| 242505 | 242505 | Rasef | chr4 | -6.539106 | 6.704195 | 2151.638 | 3.128356e-11 | 5.200579e-07 |
| 71740 | 71740 | Nectin4 | chr1 | -5.586420 | 6.447456 | 1566.715 | 1.144948e-10 | 6.188524e-07 |
| 67451 | 67451 | Pkp2 | chr16 | -5.178572 | 5.671932 | 1566.704 | 1.144981e-10 | 6.188524e-07 |
| 53624 | 53624 | Cldn7 | chr11 | -5.504748 | 7.529234 | 1333.741 | 2.210153e-10 | 6.188524e-07 |
| 78896 | 78896 | Ecrg4 | chr1 | 6.456502 | 5.672665 | 1256.914 | 2.815928e-10 | 6.188524e-07 |
| 19253 | 19253 | Ptpn18 | chr1 | -4.594745 | 5.343492 | 1155.825 | 3.964889e-10 | 6.188524e-07 |
| 50722 | 50722 | Dkkl1 | chr7 | -6.778277 | 6.097962 | 1152.189 | 4.016182e-10 | 6.188524e-07 |
| 12521 | 12521 | Cd82 | chr2 | -4.687198 | 8.007496 | 1151.233 | 4.029809e-10 | 6.188524e-07 |
| 218518 | 218518 | Marveld2 | chr13 | -4.755942 | 6.054055 | 1146.641 | 4.096058e-10 | 6.188524e-07 |
| 108153 | 108153 | Adamts7 | chr9 | 4.697586 | 5.748726 | 1135.707 | 4.259335e-10 | 6.188524e-07 |
Y ahora comparamos LP con ML.
contraste <- makeContrasts(LP - ML, levels = diseño)
glmQL_LP_ML <- glmQLFTest(modelo, contrast = contraste)
summary(decideTests(glmQL_LP_ML)) 1*LP -1*ML
Down 3175
NotSig 10879
Up 2570
as.data.frame(topTags(glmQL_LP_ML)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC | logCPM | F | PValue | FDR |
|---|---|---|---|---|---|---|---|---|
| 11815 | 11815 | Apod | chr16 | 4.283581 | 6.321010 | 486.9859 | 1.330736e-08 | 0.0002212215 |
| 20319 | 20319 | Sfrp2 | chr3 | 2.753550 | 4.598845 | 314.7591 | 7.723980e-08 | 0.0005622194 |
| 15212 | 15212 | Hexb | chr13 | 2.909540 | 5.986395 | 286.4532 | 1.126803e-07 | 0.0005622194 |
| 13132 | 13132 | Dab2 | chr15 | 2.557433 | 5.184149 | 271.5224 | 1.395871e-07 | 0.0005622194 |
| 14962 | 14962 | Cfb | chr17 | 1.803855 | 4.762677 | 253.8898 | 1.825182e-07 | 0.0005622194 |
| 18552 | 18552 | Pcsk5 | chr19 | 2.158846 | 4.236652 | 238.4773 | 2.342784e-07 | 0.0005622194 |
| 12424 | 12424 | Cck | chr9 | 4.313973 | 4.761400 | 236.2864 | 2.430493e-07 | 0.0005622194 |
| 74365 | 74365 | Lonrf3 | chrX | 2.245457 | 4.517287 | 226.4194 | 2.880104e-07 | 0.0005622194 |
| 73341 | 73341 | Arhgef6 | chrX | 2.346351 | 7.760006 | 219.2708 | 3.271778e-07 | 0.0005622194 |
| 18858 | 18858 | Pmp22 | chr11 | 1.724631 | 5.964828 | 204.3723 | 4.325480e-07 | 0.0005622194 |
Por último buscamos los genes con una expresión diferente en los tres grupos experimentales.
glmQL_all <- glmQLFTest(modelo, coef = 2:3)
as.data.frame(topTags(glmQL_all)) |>
gt(rownames_to_stub = T, caption = "Genes más subexpresados o sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | ENTREZID | SYMBOL | TXCHROM | logFC.LP | logFC.ML | logCPM | F | PValue | FDR |
|---|---|---|---|---|---|---|---|---|---|
| 68598 | 68598 | Dnajc8 | chr4 | -14.75481 | -14.52690 | 5.503444 | 3917.603 | 5.798476e-13 | 3.280446e-11 |
| 78658 | 78658 | Ncapd3 | chr9 | -15.14518 | -15.10475 | 5.096621 | 3893.545 | 5.946855e-13 | 3.280446e-11 |
| 76251 | 76251 | Ercc6l2 | chr13 | -14.63388 | -14.74224 | 5.220876 | 3886.831 | 5.989106e-13 | 3.280446e-11 |
| 67444 | 67444 | Ilkap | chr1 | -14.43347 | -14.33013 | 5.399407 | 3841.724 | 6.282843e-13 | 3.280446e-11 |
| 58859 | 58859 | Efemp2 | chr19 | -14.80616 | -14.65589 | 5.410953 | 3840.347 | 6.292087e-13 | 3.280446e-11 |
| 24045 | 24045 | Scamp3 | chr3 | -15.06093 | -14.58439 | 5.202811 | 3765.952 | 6.817772e-13 | 3.280446e-11 |
| 21371 | 21371 | Tbca | chr13 | -14.99664 | -15.05218 | 5.185991 | 3765.800 | 6.818903e-13 | 3.280446e-11 |
| 57443 | 57443 | Fbxo3 | chr2 | -14.01332 | -14.04999 | 5.861901 | 3765.551 | 6.820748e-13 | 3.280446e-11 |
| 80517 | 80517 | Herpud2 | chr9 | -14.26672 | -14.05908 | 5.853877 | 3756.947 | 6.885051e-13 | 3.280446e-11 |
| 231464 | 231464 | Cnot6l | chr5 | -14.28498 | -14.01301 | 5.835142 | 3750.520 | 6.933572e-13 | 3.280446e-11 |
1.1.5 Visualización de resultados
Existen diferentes diagramas para representar los resultados del análisis génico diferencial. En esta sección presentamos dos de los diagramas más utilizados, el diagrama MD y el diagrama de volcán.
1.1.5.1 Diagrama MD
Este diagrama consiste en un mapa de puntos donde cada punto corresponde a un gen. El eje X representa la media del logaritmo de las frecuencias por millón (lcpm) y el eje Y representa el logaritmo en base 2 de la razón de cambio (logFC). Habitualmente, los genes subexpresados significativamente se dibujan en color azul, mientras que los sobreexpresados se dibujan en color rojo.
Para dibujar un diagrama MD con edgeR se utiliza la función
plotMD(contraste): Dibuja un diagrama MD a partir de los resultados del contraste de comparación de expresión génica dado encontraste.
Ejemplo 1.14 Dibujamos primero el diagrama MD para la comparación de los grupos Basal y LP.
plotMD(glmQL_B_LP)
A continuación para la comparación de los grupos Basal y ML.
plotMD(glmQL_B_ML)
Y finalmente para la comparación de los grupos LP y ML.
plotMD(glmQL_LP_ML)
1.1.5.2 Diagrama de volcán
Un diagrama de volcán es una representación gráfica utilizada en Bioinformática para visualizar los resultados del análisis de expresión génica diferencial u otros tipos de análisis de datos ómicos de alto rendimiento, como Proteómica o Metabolómica.
En un diagrama de volcán, cada punto de datos representa un gen (o una proteína/metabolito) del conjunto de datos. El eje x muestra el cambio logarítmico o tamaño del efecto, que mide la magnitud del cambio en la expresión entre dos condiciones (por ejemplo, tratamiento vs. control). El eje y muestra la significación estadística, a menudo representada como el logaritmo negativo del valor p. Un valor p más pequeño indica una mayor significación estadística.
A menudo, los puntos en el gráfico de volcán están coloreados o resaltados según su significación estadística y cambio en la expresión. Por lo general, los genes significativamente sobreexpresados se representan en rojo, los genes significativamente subexpresados en azul y los genes no significativamente expresados en gris o negro. Cuanto más significativos sean estadísticamente y mayor sea el cambio en la expresión, más alejados estarán los puntos del centro del gráfico.
edgR no incorpora una función para realizar este tipo de diagramas, así que utilizaremos la siguiente función del paquete Glimma
glimmaVolcano(contraste, dge = dge): Realiza diagrama de volcán del contraste de comparación de expresión génica dado encontraste, realizado sobre la estructura de datos del tipoDGEListdada endge.
En el capítulo Diagramas de Volcán se explica con más detalle su creación con el paquete ggplot2.
Ejemplo 1.15 Dibujamos primero el diagrama de volcán para la comparación de los grupos Basal y LP.
library(Glimma)
glimmaVolcano(glmQL_B_LP, dge = dge)A continuación para la comparación de los grupos Basal y ML.
library(Glimma)
glimmaVolcano(glmQL_B_ML, dge = dge)Y finalmente para la comparación de los grupos LP y ML.
library(Glimma)
glimmaVolcano(glmQL_LP_ML, dge = dge)1.2 Análisis de expresión génica diferencial con DESeq2
El paquete DESeq2 es otro de los paquetes más utilizados para el análisis de la expresión génica diferencial, que al igual que edgeR incorpora procedimientos para todas las etapas del análisis. También está disponible en el repositorio Bioconductor.
1.2.1 Estructura de Datos
El paquete DESeq2 utiliza la estructura de datos DESeqDataSet para guardar la matriz de frecuencias de expresión génica. Esta estructura se deriva, a su vez, de la estructura SingleCellExperiment, que además de la matriz de frecuencias de expresión de los genes incluye también los grupos experimentales a los que pertenecen las muestras estudiadas y otra información que se va generando a medida que progresa el análisis.

SingleCellExperimentEl paquete SingleCellExperiment de Bioconductor define esta estructura de datos, y puede instalarse mediante el siguiente código
BiocManager::install('SingleCellExperiment')Para crear la estructura de datos DESeqDataSet se utiliza la siguiente función
DESeqDataSetFromMatrix(countData = frec, colData = grupo, design = diseño): Crea una estructura del tipoDESeqDataSetcon la matriz de frecuencias de expresión génicafrec(con genes en filas y muestras en columnas), y el data framegrupocon los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias y la columnadiseñodel data framegrupoque contiene los grupos experimentales a comparar.
Ejemplo 1.16 Siguiendo con el mismo ejemplo de la sección anterior, podemos aprovechar la estructura de datos DGEList creada durante el análisis génico diferencial con el paquete edgeR para crear a partir de ella la estructura de datos DESeqDataSet que requiere el paquete DESeq2.
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
countData = frecuencias,
colData = muestras,
rowData = genes,
design = ~ Grupo
)Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
design formula are characters, converting to factors
ddsclass: DESeqDataSet
dim: 27179 9
metadata(1): version
assays(1): counts
rownames(27179): 497097 100503874 ... 100861691 100504472
rowData names(3): ENTREZID SYMBOL TXCHROM
colnames(9): LP1 ML1 ... ML3 LP3
colData names(3): Id Grupo Lote
Si disponemos de los datos en el formato de la estructura de datos DGEList, es fácil convertirlos al formato de la estructura de datos DESeqDataSet mediante la función as.DESeqDataSet del paquete DEFormats.
as.DESeqDataSet(dge): Devuelve una estructura de datos del tipoDESeqDataSetcon los datos de la estructura de datosDGEListdada endge.
Ejemplo 1.17
library(DEFormats)
dds <- as.DESeqDataSet(dge)Una vez creada la estructura de datos, podemos acceder a la información que contiene con las siguientes funciones:
counts(dds): Devuelve un data frame con las frecuencias de expresión génica de la estructura de datos en formatoDESEqDataSetdada endds.colData(dds): Devuelve un data frame con la información de las muestras (columnas) de la estructura de datos en formatoDESEqDataSetdada endds.rowData(dds): Devuelve un data frame con la información de los genes (filas) de la estructura de datos en formatoDESEqDataSetdada endds.design(dds): Devuelve la fórmula del diseño experimental de la estructura de datos en formatoDESEqDataSetdada endds.
Ejemplo 1.18 Para comprobar que se creado bien la estructura de datos, mostramos la información que contiene. Primero la matriz de frecuencias de expresión génica.
head(counts(dds)) LP1 ML1 B1 B2 ML2 LP2 B3 ML3 LP3
497097 1 2 342 526 3 3 535 2 0
100503874 0 0 5 6 0 0 5 0 0
100038431 0 0 0 0 0 0 1 0 0
19888 0 1 0 0 17 2 0 1 0
20671 1 1 76 40 33 14 98 18 8
27395 431 771 1368 1268 1564 769 818 468 342
Después, la información de las muestras del ejemplo.
colData(dds)DataFrame with 9 rows and 3 columns
Id Grupo Lote
<character> <factor> <character>
LP1 LP1 LP L004
ML1 ML1 ML L004
B1 B1 Basal L004
B2 B2 Basal L006
ML2 ML2 ML L006
LP2 LP2 LP L006
B3 B3 Basal L006
ML3 ML3 ML L008
LP3 LP3 LP L008
A continuación, mostramos la información de los genes.
rowData(dds)DataFrame with 27179 rows and 3 columns
ENTREZID SYMBOL TXCHROM
<character> <character> <character>
497097 497097 Xkr4 chr1
100503874 100503874 Gm19938 NA
100038431 100038431 Gm10568 NA
19888 19888 Rp1 chr1
20671 20671 Sox17 chr1
... ... ... ...
100861837 100861837 NA NA
100861924 100861924 NA NA
170942 170942 Erdr1x chrY
100861691 100861691 LOC100861691 NA
100504472 100504472 NA NA
Y finalmente mostramos el diseño.
design(dds)~Grupo
1.2.2 Preprocesamiento
El paquete DESeq2 realiza el preprocesamiento de maneare automática cuando se lanza el procedimiento para el análisis de expresión génica diferencial.
1.2.3 Análisis exploratorio
1.2.4 Escalado multidimensional
Al igual que antes, el principal análisis exploratorio es el escalado multidimensional. Para realizarlo podemos usar de nuevo la siguiente función del paquete Glimma.
glimmaMDS(dds): Realiza un diagrama interactivo de componenetes principales a partir de la estructura de datosdds.
Ejemplo 1.19 A continuación se muestra cómo obtener el diagrama de componentes principales para nuestro ejemplo con el paquete Glimma.
library(Glimma)
glimmaMDS(dds)1.2.5 Análisis de expresión génica diferencial
Para realizar el análisis de expresión génica diferencial el paquete DESeq2 utiliza la siguiente función
DESeq(dds): Realiza el preprocesamiento de datos y el ajuste del modelo para el análisis de expresión génica diferencial de la estructura de datos en formatoDESeqDataSetdada endds.
Antes de proceder con el análisis es importante definir la categoría de referencia (grupo control) como el primer nivel del principal factor experimental. Ya que el cálculo de los estadísticos como el logaritmo en base 2 de la razón de cambio (logFC) se realizará con respecto a esta categoría.
Para visualizar los resultados de análisis de expresión génica diferencial se utilizan las funciones
results(dds, contrast = c(factor, grupo1, grupo2), alpha = 0.05): Muestra una tabla con los principales estadísticos del contraste de comparación de expresión génica delgrupo2frente algrupo1del factor experimental dado enfactor, realizado a partir de la estructura de datos en formatoDESeqDataSetdada endds. Si no se indica el parámetrocontrast, por defecto se muestra el resultado del contraste de comparación del último nivel del factor experimental frente al primero. El parámetroalphaindica el nivel de significación del contraste, des decir, probabilidad por debajo de la cuál se considera significativo el p-valor ajustado.summary(contraste): Devuelve un resumen con los genes subexpresados y sobreexpresados significativamente a partir de los resultados del contraste de comparación de expresión génica dados encontraste.
Ejemplo 1.20 Siguiendo con nuestro ejemplo, primero definiremos el grupo Basal como la primera categoría del factor group.
dds$Grupo <- relevel(dds$Grupo, "Basal")A continuación ya podemos realizar el análisis de expresión génica diferencial para nuestro ejemplo.
dds <- DESeq(dds)estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
ddsclass: DESeqDataSet
dim: 27179 9
metadata(1): version
assays(4): counts mu H cooks
rownames(27179): 497097 100503874 ... 100861691 100504472
rowData names(29): ENTREZID SYMBOL ... deviance maxCooks
colnames(9): LP1 ML1 ... ML3 LP3
colData names(4): Id Grupo Lote sizeFactor
A continuación, mostramos los resultados para la comparación del grupo Basal con el grupo LP.
res_B_LP <- results(dds, contrast = c("Grupo", "Basal", "LP"))
# Resumen de genes subexpresados y sobreexpresados.
summary(res_B_LP)
out of 22026 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 5494, 25%
LFC < 0 (down) : 5264, 24%
outliers [1] : 18, 0.082%
low counts [2] : 2953, 13%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
# Genes subexpresados significativamente ordenados de menor a mayor logFC.
as.data.frame(res_B_LP) |>
filter(log2FoldChange < 0 & padj < 0.05) |>
arrange(log2FoldChange) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más subexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 243874 | 49.26108 | -8.965084 | 1.2704518 | -7.056611 | 1.706126e-12 | 1.262043e-11 |
| 347708 | 243.18931 | -8.453980 | 1.0758631 | -7.857858 | 3.907565e-15 | 3.540592e-14 |
| 18768 | 302.95524 | -8.283820 | 0.5738121 | -14.436467 | 3.050847e-47 | 1.234265e-45 |
| 68891 | 85.08457 | -8.240424 | 0.9295045 | -8.865394 | 7.623421e-19 | 8.682862e-18 |
| 545847 | 73.96691 | -8.120701 | 0.9030947 | -8.992082 | 2.425908e-19 | 2.856963e-18 |
| 20750 | 138304.22680 | -8.046167 | 0.4702365 | -17.110896 | 1.230909e-65 | 1.037831e-63 |
| 230405 | 23.02081 | -7.939209 | 1.2914113 | -6.147700 | 7.861447e-10 | 4.715136e-09 |
| 17830 | 29.79799 | -7.915517 | 1.3274643 | -5.962885 | 2.478234e-09 | 1.412586e-08 |
| 237749 | 40.12100 | -7.847820 | 1.2804673 | -6.128871 | 8.850462e-10 | 5.283382e-09 |
| 16415 | 15.09905 | -7.834211 | 1.2996894 | -6.027757 | 1.662512e-09 | 9.673028e-09 |
# Genes sobrexpresados significativamente ordenados de mayor a menor logFC.
as.data.frame(res_B_LP) |>
filter(log2FoldChange > 0 & padj < 0.05) |>
arrange(desc(log2FoldChange)) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 624245 | 130.01315 | 10.651813 | 1.2050616 | 8.839227 | 9.638348e-19 | 1.092556e-17 |
| 360220 | 91.29897 | 9.185650 | 1.1765135 | 7.807518 | 5.832527e-15 | 5.203127e-14 |
| 100861621 | 152.72891 | 9.143254 | 1.0158933 | 9.000211 | 2.252843e-19 | 2.656431e-18 |
| 74720 | 31.74058 | 8.621050 | 1.3331640 | 6.466608 | 1.002270e-10 | 6.480572e-10 |
| 224065 | 31.39369 | 8.603568 | 1.2600963 | 6.827706 | 8.628293e-12 | 6.051237e-11 |
| 22409 | 1408.83897 | 8.488714 | 0.4741240 | 17.903996 | 1.097572e-71 | 1.130499e-69 |
| 26950 | 297.03582 | 8.157546 | 1.5000622 | 5.438138 | 5.384013e-08 | 2.687071e-07 |
| 100861755 | 228.37922 | 7.958092 | 0.6120430 | 13.002505 | 1.184006e-38 | 3.444462e-37 |
| 23964 | 4653.79258 | 7.941830 | 0.4989370 | 15.917500 | 4.791431e-57 | 2.880149e-55 |
| 223838 | 501.88062 | 7.843173 | 0.4330005 | 18.113540 | 2.491951e-73 | 2.728972e-71 |
Después, mostramos los resultados para la comparación del grupo Basal con el grupo ML.
res_B_ML <- results(dds, contrast = c("Grupo", "Basal", "ML"))
# Resumen de genes subexpresados y sobreexpresados.
summary(res_B_ML)
out of 22026 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 5618, 26%
LFC < 0 (down) : 5402, 25%
outliers [1] : 18, 0.082%
low counts [2] : 2953, 13%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
# Genes subexpresados significativamente ordenados de menor a mayor logFC.
as.data.frame(res_B_LP) |>
filter(log2FoldChange < 0 & padj < 0.05) |>
arrange(log2FoldChange) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más subexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 243874 | 49.26108 | -8.965084 | 1.2704518 | -7.056611 | 1.706126e-12 | 1.262043e-11 |
| 347708 | 243.18931 | -8.453980 | 1.0758631 | -7.857858 | 3.907565e-15 | 3.540592e-14 |
| 18768 | 302.95524 | -8.283820 | 0.5738121 | -14.436467 | 3.050847e-47 | 1.234265e-45 |
| 68891 | 85.08457 | -8.240424 | 0.9295045 | -8.865394 | 7.623421e-19 | 8.682862e-18 |
| 545847 | 73.96691 | -8.120701 | 0.9030947 | -8.992082 | 2.425908e-19 | 2.856963e-18 |
| 20750 | 138304.22680 | -8.046167 | 0.4702365 | -17.110896 | 1.230909e-65 | 1.037831e-63 |
| 230405 | 23.02081 | -7.939209 | 1.2914113 | -6.147700 | 7.861447e-10 | 4.715136e-09 |
| 17830 | 29.79799 | -7.915517 | 1.3274643 | -5.962885 | 2.478234e-09 | 1.412586e-08 |
| 237749 | 40.12100 | -7.847820 | 1.2804673 | -6.128871 | 8.850462e-10 | 5.283382e-09 |
| 16415 | 15.09905 | -7.834211 | 1.2996894 | -6.027757 | 1.662512e-09 | 9.673028e-09 |
# Genes sobrexpresados significativamente ordenados de mayor a menor logFC.
as.data.frame(res_B_LP) |>
filter(log2FoldChange > 0 & padj < 0.05) |>
arrange(desc(log2FoldChange)) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 624245 | 130.01315 | 10.651813 | 1.2050616 | 8.839227 | 9.638348e-19 | 1.092556e-17 |
| 360220 | 91.29897 | 9.185650 | 1.1765135 | 7.807518 | 5.832527e-15 | 5.203127e-14 |
| 100861621 | 152.72891 | 9.143254 | 1.0158933 | 9.000211 | 2.252843e-19 | 2.656431e-18 |
| 74720 | 31.74058 | 8.621050 | 1.3331640 | 6.466608 | 1.002270e-10 | 6.480572e-10 |
| 224065 | 31.39369 | 8.603568 | 1.2600963 | 6.827706 | 8.628293e-12 | 6.051237e-11 |
| 22409 | 1408.83897 | 8.488714 | 0.4741240 | 17.903996 | 1.097572e-71 | 1.130499e-69 |
| 26950 | 297.03582 | 8.157546 | 1.5000622 | 5.438138 | 5.384013e-08 | 2.687071e-07 |
| 100861755 | 228.37922 | 7.958092 | 0.6120430 | 13.002505 | 1.184006e-38 | 3.444462e-37 |
| 23964 | 4653.79258 | 7.941830 | 0.4989370 | 15.917500 | 4.791431e-57 | 2.880149e-55 |
| 223838 | 501.88062 | 7.843173 | 0.4330005 | 18.113540 | 2.491951e-73 | 2.728972e-71 |
Y finalmente mostramos los resultados para la comparación del grupo LP con el grupo ML.
res_LP_ML <- results(dds, contrast = c("Grupo", "LP", "ML"))
# Resumen de genes subexpresados y sobreexpresados.
summary(res_LP_ML)
out of 22026 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 3268, 15%
LFC < 0 (down) : 3471, 16%
outliers [1] : 18, 0.082%
low counts [2] : 3796, 17%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
# Genes subexpresados significativamente ordenados de menor a mayor logFC.
as.data.frame(res_LP_ML) |>
filter(log2FoldChange < 0 & padj < 0.05) |>
arrange(log2FoldChange) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más subexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 16675 | 35.256494 | -6.365421 | 2.054093 | -3.098895 | 1.942436e-03 | 0.0092244161 |
| 70843 | 73.553279 | -6.295554 | 1.593891 | -3.949803 | 7.821551e-05 | 0.0006064116 |
| 116746 | 10.174747 | -5.880935 | 1.945496 | -3.022846 | 2.504096e-03 | 0.0113813883 |
| 20423 | 4.821526 | -5.569531 | 1.664008 | -3.347058 | 8.167400e-04 | 0.0045238652 |
| 13392 | 38.113342 | -5.423817 | 1.826760 | -2.969090 | 2.986827e-03 | 0.0130980234 |
| 22415 | 7.388444 | -5.378074 | 1.853749 | -2.901187 | 3.717518e-03 | 0.0156975269 |
| 319180 | 6.635441 | -5.256345 | 1.599572 | -3.286094 | 1.015873e-03 | 0.0054160042 |
| 19293 | 5.787450 | -5.185556 | 1.477070 | -3.510703 | 4.469229e-04 | 0.0027095072 |
| 13190 | 6.282620 | -5.102299 | 2.006739 | -2.542582 | 1.100367e-02 | 0.0382542306 |
| 242627 | 5.890588 | -5.038786 | 1.904499 | -2.645727 | 8.151545e-03 | 0.0299185702 |
# Genes sobrexpresados significativamente ordenados de mayor a menor logFC.
as.data.frame(res_LP_ML) |>
filter(log2FoldChange > 0 & padj < 0.05) |>
arrange(desc(log2FoldChange)) |>
head(10) |>
gt(rownames_to_stub = T, caption = "10 genes más sobreexpresados signiticativamente") |>
tab_stubhead(label = "Gen")| Gen | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| 81906 | 5.084982e+00 | 6.349445 | 1.5596581 | 4.071049 | 4.680186e-05 | 3.900940e-04 |
| 83563 | 8.780151e+00 | 5.599738 | 1.7943674 | 3.120731 | 1.804026e-03 | 8.696378e-03 |
| 75556 | 2.880224e+00 | 5.504546 | 1.8676200 | 2.947359 | 3.205013e-03 | 1.389755e-02 |
| 67719 | 8.400270e+00 | 4.993597 | 1.2677240 | 3.939025 | 8.181335e-05 | 6.297484e-04 |
| 347708 | 2.431893e+02 | 4.966281 | 0.9673189 | 5.134068 | 2.835461e-07 | 4.451674e-06 |
| 245446 | 1.262976e+02 | 4.750196 | 0.4521869 | 10.504939 | 8.197556e-26 | 1.492939e-22 |
| 100470 | 2.577426e+03 | 4.594474 | 1.1793063 | 3.895912 | 9.782983e-05 | 7.365345e-04 |
| 12993 | 3.527468e+04 | 4.581113 | 1.1971312 | 3.826742 | 1.298503e-04 | 9.373102e-04 |
| 20750 | 1.383042e+05 | 4.541280 | 0.4699195 | 9.663952 | 4.289877e-22 | 2.441476e-19 |
| 14395 | 3.024909e+00 | 4.528923 | 1.5358977 | 2.948714 | 3.190989e-03 | 1.384663e-02 |
1.2.6 Visualización de datos
1.2.6.1 Diagrama MD
Este diagrama consiste en un mapa de puntos donde cada punto corresponde a un gen. El eje X representa la media del logaritmo de las frecuencias por millón (lcpm) y el eje Y representa el logaritmo en base 2 de la razón de cambio (logFC). Habitualmente, los genes subexpresados significativamente se dibujan en color azul, mientras que los sobreexpresados se dibujan en color rojo.
DESeq2 no incorpora una función para realizar este tipo de diagramas, así que utilizaremos la siguiente función del paquete Glimma
glimmaMA(dds, groups = factor): Realiza diagrama MD del contraste de comparación de expresión génica realizado sobre la estructura de datos del tipoDESeqDataSetdada endds, con los grupos experimentales definidos enfactor.
La función glimmaMA muestra el diagrama de comparación para la comparación del último nivel del factor experimental con el primero. Si tenemos más de dos grupos experimentales hay que filtrar previamente la estructura de datos en formato DESeqDataSet para que solo incluya los dos grupos a comparar.
Ejemplo 1.21
glimmaMA(dds, groups = dds$Grupo)8124 of 27179 genes were filtered out in DESeq2 tests
1.2.6.2 Diagrama de volcán
Un diagrama de volcán es una representación gráfica utilizada en Bioinformática para visualizar los resultados del análisis de expresión génica diferencial u otros tipos de análisis de datos ómicos de alto rendimiento, como Proteómica o Metabolómica.
En un diagrama de volcán, cada punto de datos representa un gen (o una proteína/metabolito) del conjunto de datos. El eje x muestra el cambio logarítmico o tamaño del efecto, que mide la magnitud del cambio en la expresión entre dos condiciones (por ejemplo, tratamiento vs. control). El eje y muestra la significación estadística, a menudo representada como el logaritmo negativo del valor p. Un valor p más pequeño indica una mayor significación estadística.
A menudo, los puntos en el gráfico de volcán están coloreados o resaltados según su significación estadística y cambio en la expresión. Por lo general, los genes significativamente sobreexpresados se representan en rojo, los genes significativamente subexpresados en azul y los genes no significativamente expresados en gris o negro. Cuanto más significativos sean estadísticamente y mayor sea el cambio en la expresión, más alejados estarán los puntos del centro del gráfico.
DESeq2 no incorpora una función para realizar este tipo de diagramas, pero podemos usar la siguiente función del paquete Glimma
glimmaVolcano(dds, groups = factor): Realiza diagrama de volcán del contraste de comparación de expresión génica realizado sobre la estructura de datos del tipoDESeqDataSetdada endds, con los grupos experimentales definidos enfactor.
La función glimmaVolcano muestra el diagrama de comparación para la comparación del último nivel del factor experimental con el primero. Si tenemos más de dos grupos experimentales hay que filtrar previamente la estructura de datos en formato DESeqDataSet para que solo incluya los dos grupos a comparar.
En el capítulo Diagramas de Volcán se explica con más detalle su creación con el paquete ggplot2.
Ejemplo 1.22 Dibujamos el diagrama de volcán interactivo.
glimmaVolcano(dds, groups = dds$Grupo)Obsérvese que este diagrama solo muestra la comparación del último grupo experimental (ML) con respecto al grupo de referencia (Basal).
A continuación creamos los diagrams de volcán para todos los contrastes por pares, mediante el paquete ggplot. En primer lugar para la comparación de LP frente a Basal.
# Crear nueva columna con categorías de expresión génica.
diagramaVolcan <- function(contraste, alpha = 0.05, minlogFC = 2){
colores <- c("dodgerblue3", "gray50", "firebrick3")
df <- as.data.frame(contraste) |> na.omit()
rangex <- max(-min(df$log2FoldChange, na.rm = T), max(df$log2FoldChange, na.rm = T))
maxy <- max(-log10(df$padj), na.rm = T)
vplot <- df |>
mutate(Expresión = case_when(
log2FoldChange >= log2(minlogFC) & padj <= alpha ~ "Up",
log2FoldChange <= log2(1/minlogFC) & padj <= alpha ~ "Down",
TRUE ~ "NC")) |>
# Definimos los colores para los genes subexpresados, normales y sobreexpresados.
ggplot(aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = Expresión)) +
scale_color_manual(values = colores) +
# Añadimos las líneas de la razón de cambio.
geom_vline(xintercept = c(log2(1/minlogFC), log2(minlogFC)), linetype = "dashed") +
# Etiquetamos las líneas de la razón de cambio.
annotate("text", x = log2(1/minlogFC), y = maxy, label = paste0("FC=1/", minlogFC), vjust = 0) +
annotate("text", x = log2(2), y = maxy, label = paste0("FC=", minlogFC), vjust = 0) +
# Añadimos las líneas de los p-valores.
geom_hline(yintercept = -log10(alpha), linetype = "dashed") +
# Etiquetamos las líneas de los p-valores.
annotate("text", x = rangex+0.1, y = -log10(alpha), label = paste0("p=", alpha), hjust = 0) +
theme(plot.margin = unit(c(2, 1, 1, 1), "lines")) +
coord_cartesian(xlim = c(-rangex, rangex), ylim = c(0, maxy), expand = FALSE, clip = "off")
return(vplot)
}diagramaVolcan(res_B_LP, alpha = 0.05, minlogFC = 2)
A continuación dibujamos el diagrama de volcán de ML frente a Basal.
diagramaVolcan(res_B_ML, alpha = 0.05, minlogFC = 2)
Y finalmente el diagrama de volcán de ML frente a LP.
diagramaVolcan(res_LP_ML, alpha = 0.05, minlogFC = 2)
1.3 Referencias
Orchestrating Single-Cell Analysis with Bioconductor. (s. f.). Recuperado 12 de junio de 2023, de https://github.com/LTLA/OSCA
Law, Charity, et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR.